home *** CD-ROM | disk | FTP | other *** search
/ PC Media 23 / PC MEDIA CD23.iso / share / prog / newmat / newmat5.cpp < prev    next >
C/C++ Source or Header  |  1995-01-17  |  12KB  |  435 lines
  1. //$$ newmat5.cpp         Transpose, evaluate etc
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,5); ++ExeCount; }
  11.  
  12. #define REPORT {}
  13.  
  14.  
  15. /************************ carry out operations ******************************/
  16.  
  17.  
  18. GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
  19. {
  20.    GeneralMatrix* gm1;
  21.  
  22.    if (Compare(Type().t(),mt))
  23.    {
  24.       REPORT
  25.       gm1 = mt.New(ncols,nrows,tm);
  26.       for (int i=0; i<ncols; i++)
  27.       {
  28.      MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
  29.          MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
  30.       }
  31.    }
  32.    else
  33.    {
  34.       REPORT
  35.       gm1 = mt.New(ncols,nrows,tm);
  36.       MatrixRow mr(this, LoadOnEntry);
  37.       MatrixCol mc(gm1, StoreOnExit+DirectPart);
  38.       int i = nrows;
  39.       while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
  40.    }
  41.    tDelete(); gm1->ReleaseAndDelete(); return gm1;
  42. }
  43.  
  44. GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  45. { REPORT  return Evaluate(mt); }
  46.  
  47.  
  48. GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  49. { REPORT return Evaluate(mt); }
  50.  
  51. Boolean GeneralMatrix::IsZero() const
  52. {
  53.    REPORT
  54.    Real* s=store; int i=storage;
  55.    while (i--) { if (*s++) return FALSE; }
  56.    return TRUE;
  57. }
  58.  
  59. GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
  60. {
  61.    REPORT
  62.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  63.    gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
  64.    return BorrowStore(gmx,mt);
  65. }
  66.  
  67. GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
  68. {
  69.    REPORT
  70.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  71.    gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
  72.    return BorrowStore(gmx,mt);
  73. }
  74.  
  75. GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
  76. {
  77.    if (Compare(this->Type(),mt)) { REPORT return this; }
  78.    REPORT
  79.    GeneralMatrix* gmx = mt.New(nrows,ncols,this);
  80.    MatrixRow mr(this, LoadOnEntry);
  81.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  82.    int i=nrows;
  83.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  84.    tDelete(); gmx->ReleaseAndDelete(); return gmx;
  85. }
  86.  
  87. GeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
  88. {
  89.    if (Compare(cgm->Type(),mt))
  90.    {
  91.       REPORT
  92. #ifdef TEMPS_DESTROYED_QUICKLY
  93.       GeneralMatrix* gmx = (GeneralMatrix*)cgm; delete this; return gmx;
  94. #else
  95.       return (GeneralMatrix*)cgm;
  96. #endif
  97.    }
  98.    REPORT
  99.    GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols(),this);
  100.    MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
  101.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  102.    int i=cgm->Nrows();
  103.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  104.    gmx->ReleaseAndDelete();
  105. #ifdef TEMPS_DESTROYED_QUICKLY
  106.    delete this;
  107. #endif
  108.    return gmx; // no tDelete
  109. }
  110.  
  111. GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
  112. {
  113.    gm=((BaseMatrix*&)bm)->Evaluate();
  114.    int nr=gm->Nrows(); int nc=gm->Ncols();
  115.    Compare(gm->Type().AddEqualEl(),mt);
  116.    if (!(mt==gm->Type()))
  117.    {
  118.       REPORT
  119.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  120.       MatrixRow mr(gm, LoadOnEntry); 
  121.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  122.       while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
  123.       gmx->ReleaseAndDelete(); gm->tDelete();
  124. #ifdef TEMPS_DESTROYED_QUICKLY
  125.       delete this;
  126. #endif
  127.       return gmx;
  128.    }
  129.    else if (gm->reuse())
  130.    {
  131.       REPORT gm->Add(f);
  132. #ifdef TEMPS_DESTROYED_QUICKLY
  133.       GeneralMatrix* gmx = gm; delete this; return gmx;
  134. #else
  135.       return gm;
  136. #endif
  137.    }
  138.    else
  139.    {
  140.       REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
  141.       gmy->ReleaseAndDelete(); gmy->Add(gm,f);
  142. #ifdef TEMPS_DESTROYED_QUICKLY
  143.       delete this;
  144. #endif
  145.       return gmy;
  146.    }
  147. }
  148.  
  149. GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
  150. {
  151.    gm=((BaseMatrix*&)bm)->Evaluate();
  152.    int nr=gm->Nrows(); int nc=gm->Ncols();
  153.    if (Compare(gm->Type(),mt))
  154.    {
  155.       if (gm->reuse())
  156.       {
  157.          REPORT gm->Multiply(f);
  158. #ifdef TEMPS_DESTROYED_QUICKLY
  159.          GeneralMatrix* gmx = gm; delete this; return gmx;
  160. #else
  161.          return gm;
  162. #endif
  163.       }
  164.       else
  165.       {
  166.          REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  167.          gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
  168. #ifdef TEMPS_DESTROYED_QUICKLY
  169.          delete this;
  170. #endif
  171.          return gmx;
  172.       }
  173.    }
  174.    else
  175.    {
  176.       REPORT
  177.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  178.       MatrixRow mr(gm, LoadOnEntry); 
  179.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  180.       while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
  181.       gmx->ReleaseAndDelete(); gm->tDelete();
  182. #ifdef TEMPS_DESTROYED_QUICKLY
  183.       delete this;
  184. #endif
  185.       return gmx;
  186.    }
  187. }
  188.  
  189. GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
  190. {
  191.    gm=((BaseMatrix*&)bm)->Evaluate();
  192.    int nr=gm->Nrows(); int nc=gm->Ncols();
  193.    if (Compare(gm->Type(),mt))
  194.    {
  195.       if (gm->reuse())
  196.       {
  197.          REPORT gm->Negate();
  198. #ifdef TEMPS_DESTROYED_QUICKLY
  199.          GeneralMatrix* gmx = gm; delete this; return gmx;
  200. #else
  201.          return gm;
  202. #endif
  203.       }
  204.       else
  205.       {
  206.          REPORT
  207.          GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  208.          gmx->ReleaseAndDelete(); gmx->Negate(gm);
  209. #ifdef TEMPS_DESTROYED_QUICKLY
  210.          delete this;
  211. #endif
  212.          return gmx;
  213.       }
  214.    }
  215.    else
  216.    {
  217.       REPORT
  218.       GeneralMatrix* gmx = mt.New(nr,nc,this);
  219.       MatrixRow mr(gm, LoadOnEntry); 
  220.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  221.       while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
  222.       gmx->ReleaseAndDelete(); gm->tDelete();
  223. #ifdef TEMPS_DESTROYED_QUICKLY
  224.       delete this;
  225. #endif
  226.       return gmx;
  227.    }
  228. }   
  229.  
  230. GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
  231. {
  232.    REPORT
  233.    gm=((BaseMatrix*&)bm)->Evaluate();
  234.    Compare(gm->Type().t(),mt);
  235.    GeneralMatrix* gmx=gm->Transpose(this, mt);
  236. #ifdef TEMPS_DESTROYED_QUICKLY
  237.    delete this;
  238. #endif
  239.    return gmx;
  240. }
  241.    
  242. GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
  243. {
  244.    gm = ((BaseMatrix*&)bm)->Evaluate();
  245.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  246.    gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
  247. #ifdef TEMPS_DESTROYED_QUICKLY
  248.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  249. #else
  250.    return gm->BorrowStore(gmx,mt);
  251. #endif
  252. }
  253.  
  254. GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
  255. {
  256.    gm = ((BaseMatrix*&)bm)->Evaluate();
  257.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  258.    gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
  259. #ifdef TEMPS_DESTROYED_QUICKLY
  260.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  261. #else
  262.    return gm->BorrowStore(gmx,mt);
  263. #endif
  264. }
  265.  
  266. GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
  267. {
  268.    gm = ((BaseMatrix*&)bm)->Evaluate();
  269.    GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
  270.    gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
  271. #ifdef TEMPS_DESTROYED_QUICKLY
  272.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  273. #else
  274.    return gm->BorrowStore(gmx,mt);
  275. #endif
  276. }
  277.  
  278. GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
  279. {
  280.    Tracer tr("MatedMatrix::Evaluate");
  281.    gm = ((BaseMatrix*&)bm)->Evaluate();
  282.    GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
  283.    gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
  284.    if (nr*nc != gmx->storage)
  285.       Throw(IncompatibleDimensionsException());
  286. #ifdef TEMPS_DESTROYED_QUICKLY
  287.    GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  288. #else
  289.    return gm->BorrowStore(gmx,mt);
  290. #endif
  291. }
  292.  
  293. GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
  294. {
  295.    REPORT
  296.    Tracer tr("SubMatrix(evaluate)");
  297.    gm = ((BaseMatrix*&)bm)->Evaluate();
  298.    if (row_number < 0) row_number = gm->Nrows();
  299.    if (col_number < 0) col_number = gm->Ncols();
  300.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  301.       Throw(SubMatrixDimensionException());
  302.    if (IsSym) Compare(gm->Type().ssub(), mt); 
  303.    else Compare(gm->Type().sub(), mt);
  304.    GeneralMatrix* gmx = mt.New(row_number, col_number, this);
  305.    int i = row_number;
  306.    MatrixRow mr(gm, LoadOnEntry, row_skip); 
  307.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  308.    MatrixRowCol sub;
  309.    while (i--)
  310.    {
  311.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  312.       mrx.Copy(sub); mrx.Next(); mr.Next();
  313.    }
  314.    gmx->ReleaseAndDelete(); gm->tDelete();
  315. #ifdef TEMPS_DESTROYED_QUICKLY
  316.    delete this;
  317. #endif
  318.    return gmx;
  319. }   
  320.  
  321.  
  322. GeneralMatrix* ReturnMatrixX::Evaluate(MatrixType mt)
  323. {
  324. #ifdef TEMPS_DESTROYED_QUICKLY_R
  325.    GeneralMatrix* gmx = gm; delete this; return gmx->Evaluate(mt);
  326. #else
  327.    return gm->Evaluate(mt);
  328. #endif
  329. }
  330.  
  331.  
  332. void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
  333. {
  334.    REPORT
  335.    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  336.    while (i--)
  337.    { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
  338.    i = storage & 3; while (i--) *s++ = *s1++ + f;
  339. }
  340.    
  341. void GeneralMatrix::Add(Real f)
  342. {
  343.    REPORT
  344.    Real* s=store; int i=(storage >> 2);
  345.    while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
  346.    i = storage & 3; while (i--) *s++ += f;
  347. }
  348.    
  349. void GeneralMatrix::Negate(GeneralMatrix* gm1)
  350. {
  351.    // change sign of elements
  352.    REPORT
  353.    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  354.    while (i--)
  355.    { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
  356.    i = storage & 3; while(i--) *s++ = -(*s1++);
  357. }
  358.    
  359. void GeneralMatrix::Negate()
  360. {
  361.    REPORT
  362.    Real* s=store; int i=(storage >> 2);
  363.    while (i--)
  364.    { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
  365.    i = storage & 3; while(i--) { *s = -(*s); s++; }
  366. }
  367.    
  368. void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
  369. {
  370.    REPORT
  371.    Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
  372.    while (i--)
  373.    { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
  374.    i = storage & 3; while (i--) *s++ = *s1++ * f;
  375. }
  376.    
  377. void GeneralMatrix::Multiply(Real f)
  378. {
  379.    REPORT
  380.    Real* s=store; int i=(storage >> 2);
  381.    while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
  382.    i = storage & 3; while (i--) *s++ *= f;
  383. }
  384.    
  385.  
  386. /************************ MatrixInput routines ****************************/
  387.  
  388. int MatrixInput::n;           // number values still to be read
  389. Real* MatrixInput::r;         // pointer to last thing read
  390. int MatrixInput::depth;       // number of objects of this class in existence
  391.  
  392. MatrixInput MatrixInput::operator<<(Real f)
  393. {
  394.    if (!(n--))
  395.    { depth=0;  n=0; Throw(ProgramException("List of values too long")); }
  396.    *(++r) = f;
  397.    return MatrixInput();
  398. }
  399.  
  400.  
  401. MatrixInput BandMatrix::operator<<(Real)
  402. {
  403.    Throw(ProgramException("Cannot use list read with a BandMatrix"));
  404.    return MatrixInput();
  405. }
  406.  
  407. void BandMatrix::operator<<(const Real*)
  408. { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
  409.  
  410. MatrixInput GeneralMatrix::operator<<(Real f)
  411. {
  412.    if (MatrixInput::n)
  413.    {
  414.       MatrixInput::depth=0;            // so we can recover
  415.       MatrixInput::n=0;                // so we can recover
  416.       Throw(ProgramException("A list of values was too short"));
  417.    }
  418.    MatrixInput::n = Storage();
  419.    if (MatrixInput::n<=0)
  420.       Throw(ProgramException("Loading data to zero dimension matrix"));
  421.    MatrixInput::r = Store(); *(MatrixInput::r) = f; MatrixInput::n--;
  422.    return MatrixInput(); 
  423. }
  424.  
  425. MatrixInput::~MatrixInput()
  426. {
  427.    if (depth  == 1 && n != 0)
  428.    {
  429.       depth = 0; n = 0;                // so we can recover
  430.       Throw(ProgramException("A list of values was too short"));
  431.     }
  432.     else if (depth>0) depth--;
  433. }
  434.  
  435.